clear;
exec('circshift.sci');
stacksize('max');

N = 500;
Na = 4;
fs = 43e6;
ts = 1/fs;
fddc = 7e6;
//%%%%%%%%%%%%%%%%% формирование сигналов %%%%%%%%%%%%%%%%%

 Xadc = csvRead('D:\Program Files\ADC_data\Agl_1601_5_minus33.csv');
//Xadc = csvRead('Xadc20.csv');
//[x1, x2, x3, x4] = read_hex_data('D:\Program Files\ADC_data\nois3_60dBm_GN22\',1);
//Xadc = [x1; x2; x3; x4];
[num, len] = size(Xadc(:,2:$));

x = Xadc(2:2:$,2:$);
//x = [Xadc(1:2:$,2:$); Xadc(2:2:$,2:$)];// I,I,I,I,Q,Q,Q,Q

t = ts:ts:len*ts;
i = sin(2*%pi*fddc*t);
q = cos(2*%pi*fddc*t);
order = 8;
[wft, wfm, fr] = wfir('lp',order,[0.25 0], 're',[3 1.8]);
//wft = hilb(25);

xi = [];
xq = [];
for n=1:Na
    m = mean(x(n,:));
    
//    xi = [xi; x(n,:) - mean(x(n,:))];
//    xq = [xq; x(n+Na,:) - mean(x(n+Na,:))];

    xi = [xi; filter(wft, 1, (x(n,:) - m).*i)];
    xq = [xq; filter(wft, 1, (x(n,:) - m).*q)];
//    
//    xi = [xi; x(n,:) - m];
//    xq = [xq; filter(wft, 1, x(n,:) - m)];
  
//    xh = hilbert(x(n,:)-m);
//    xi = [xi; real(xh)];
//    xq = [xq; imag(xh)];
end

y = [xi; xq];

//%%%%%%%%%%%%%%%%%%%% MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
main_ch = 2;
main = [main_ch main_ch+Na];
E = y(main,:);
//%  ???? ????????
//% Z = [x ; circshift(x',1)'; circshift(x',2)'; circshift(x',3)'; circshift(x',4)'; circshift(x',5)'];
//% Z(main,:) = [];
Z_data = y;
Z_data(main,:) = [];
n_delay = 3;
Z = [];

for i = 0:n_delay
    Z = [Z; circshift(Z_data',i)'];
end

//Z_ = Z_data(:,1:N);
//
//for i = 1:2:n_delay
//    Z_ = [Z_ ; Z( (i*(num-2)+1) : (i+1)*(num-2), 1:N)];
//    if i == n_delay then
//        break;
//    end
//    tmp = circshift(Z_data(:,1:N)',(i+1))';
//    Z_ = [Z_ ; tmp];
//end

db = [];
stdEout = [];
db_ = [];
stdEout_ = [];
wm = zeros(2,2*(Na - 1)*(n_delay + 1));
wm_ = zeros(2,2*(Na - 1)*(n_delay + 1));
step_iter = 100;
iter = 1;
for i = 1:step_iter:size(y,'c')-N
    start = i;//grand(1,1,'uin',1,size(y,'c')-N);
    T = start:start+N-1
    E_ = E(:,i:N+i-1);
    Z_ = Z(:,i:N+i-1);
    out = E_ - wm*Z_;

    Rxx = Z_*Z_';
    Rxd = out*Z_';
    //pause
    invRxx = inv(Rxx);
    
    wm = wm + Rxd*invRxx;
    Eout = E(:,T) - wm*Z(:,T);
    
    stdEout = [stdEout, stdev(Eout, 'c')];
    db = [db, 20*log10(stdev(E_, 'c')./stdEout(:,iter))];

//    Rxx = Z_*Z_';
//    Rxd = E_*Z_';
//    //pause
//    invRxx = inv(Rxx);
//    wm_ = Rxd*invRxx;
//    Eout_ = E(:,T) - wm_*Z(:,T);
//    
//    stdEout_ = [stdEout_, stdev(Eout_, 'c')];
//    db_ = [db_, 20*log10(stdev(E_, 'c')./stdEout_(:,iter))];
    
    iter = iter + 1;

end

figure();
plot(db');
//figure();
//plot(db_');
